Generic Chemistry
Generic chemical equilibrium and kinetic reactions following the standard approach used by reaction-transport codes such as PHREEQ and CrunchFlow, see eg (Steefel et al., 2015). This exploits timescale separation between "fast" (assumed instantaneous) chemical equilibrium reactions, and "slow" kinetic reactions or transport.
- The chemical system is represented by a small number of
totals
(or components
) and an equal number of primary
species concentrations, with secondary
species concentrations calculated from the primary species via a set of equilibrium reactions. Primary species concentrations are determined by solving the set of algebraic equations given by the constraints on total concentrations.
- Kinetic reactions (with any species,
primary
orsecondary
as reactants) are then written withtotals
as products. - Bulk transport (eg ocean advection or eddy diffusivity) transports
totals
. Molecular diffusivity (eg in a sediment) transportsprimary
orsecondary
species and accumulates fluxes intototals
.
Reservoirs
PALEOaqchem.ReservoirsAq.ReactionConstraintReservoir
— TypeReactionConstraintReservoir
A primary species and (algebraic) constraint on a corresponding total or component.
The primary species concentration or amount is defined as a PALEO State
Variable, which depending on the primary_variable
parameter, may be:
Primary_conc
: (mol m-3)Primary
: (mol)Primary_pconc
: -log 10 (concentration (mol kg-1))
The corresponding R_constraint_conc
or R_constraint
(mol) defining the algebraic constraint on the corresponding total (for use by the numerical solver) is defined as a PALEO Constraint
Variable.
This ReactionConstraintReservoir
would usually be used in combination with a ReactionReservoir
that provides the required total component concentration or amount as an ODE variable (where as usual reaction source and sink fluxes are applied to the corresponding _sms
variable). Depending on the constraint_variable
parameter, the total component may be supplied as either a per-cell concentration or amount:
R_conc
: (mol m-3)R
: (mol)
Equilibrium reactions defining secondary species should add their contributions to the total to R_calc
(mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume
is added to R_calc
(where for the usual case parameter primary_total_stoich
should be set to 1.0
). Primary species contributions to other totals can be included by setting the primary_other_components
parameter.
The numerical solver then solves for the primary species (and hence the secondary species concentrations) that (depending on the constraint_variable
parameter) satisfy one of:
0 = R_constraint_conc = R_conc - R_calc/volume
0 = R_constraint = R - R_calc
Volume conversions
The total species concentration R_conc
and primary species concentration Primary_conc
use (potentially different) volume conversions provided in volume
and primary_volume
respectively.
This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc
refers to a cell total volume, and Primary_conc
to a solute concentration.
Parameters
primary_total_stoich[Float64]
=1.0,default_value
=1.0,description
="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"primary_other_components[Vector{String}]
=String[],default_value
=String[],description
="contribution of primary species to other element or component total concentrations"primary_variable[String]
="concentration",default_value
="concentration",allowed_values
=["concentration", "amount", "p_concentration"],description
="units for primary variable"constraint_variable[String]
="concentration",default_value
="concentration",allowed_values
=["concentration", "amount"],description
="units for constraint variable"
Methods and Variables for default Parameters
do_constraintreservoir_primary
R_calc
(mol),VT_ReactContributor
,description
="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"primary_volume
–> volume (m3),VT_ReactDependency
,description
="cell volume (as used by Primary_conc)"Primary_conc
(mol m-3),VT_ReactDependency
,VF_State
,description
="concentration of primary species"
do_constraintreservoir_constraint
R_calc
(mol),VT_ReactTarget
,description
="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"R_constraint_conc
(mol m-3),VT_ReactContributor
,VF_Constraint
,description
="algebraic constraint on R_conc (= 0)"R_conc
(mol m-3),VT_ReactDependency
,description
="total R_conc"volume
(m3),VT_ReactDependency
,description
="cell volume (as used by total variable)"
PALEOaqchem.ReservoirsAq.ReactionImplicitReservoir
— TypeReactionImplicitReservoir
A primary species and corresponding total or component as an 'implicit' ODE variable.
This provides an implementation of the 'Direct Substitution Approach' to chemical speciation, where the total or component is a function of the primary species concentration.
The primary species concentration or amount is defined as a PALEO StateTotal
Variable, which depending on the primary_variable
parameter, may be:
Primary_conc
: (mol m-3)Primary
: (mol)Primary_pconc
: -log 10 (concentration (mol kg-1))
The corresponding total component R_conc
or R
is defined as a PALEO Total
Variable, which depending on the constraint_variable
parameter, may be provided to the solver either as a per-cell concentration or amount:
R_conc = R_calc/volume
: (mol m-3)R = R_calc
: (mol)
Equilibrium reactions defining secondary species should add their contributions to the total to R_calc
(mol). A primary species contribution R_calc += primary_total_stoich * Primary_conc * primary_volume
is added to R_calc
(where for the usual case parameter primary_total_stoich
should be set to 1.0
). Primary species contributions to other totals can be included by setting the primary_other_components
parameter.
Source - sink fluxes eg kinetic reactions should be added to R_sms
(mol yr-1) defined as a PALEO Deriv
Variable.
Volume conversions
The total species concentration R_conc
and primary species concentration Primary_conc
use (potentially different) volume conversions provided in volume
and primary_volume
respectively.
This allows for cases eg equilibrium partitioning between solute and solid phases by surface complexation, where R_conc
refers to a cell total volume, and Primary_conc
to a solute concentration.
Parameters
primary_total_stoich[Float64]
=1.0,default_value
=1.0,description
="stoichiometric factor R_calc_conc += primary_total_stoich * Primary_conc"primary_other_components[Vector{String}]
=String[],default_value
=String[],description
="contribution of primary species to other element or component total concentrations"primary_variable[String]
="concentration",default_value
="concentration",allowed_values
=["concentration", "amount", "pconcentration"],description
="units for primary variable (specifies Primary\conc, Primary, Primary_pconc as StateTotal variable)"total_variable[String]
="concentration",default_value
="concentration",allowed_values
=["concentration", "amount"],description
="units for total variable (specifies R_conc, R as Total variable)"total[Bool]
=false,default_value
=false,description
="true to calculate R_total = sum(R)"
Methods and Variables for default Parameters
do_implicitreservoir_primary
R_calc
(mol),VT_ReactContributor
,description
="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"primary_volume
–> volume (m3),VT_ReactDependency
,description
="cell volume (as used by Primary_conc)"Primary_conc
(mol m-3),VT_ReactDependency
,VF_StateTotal
,description
="concentration of primary species"
do_implicitreservoir_sms
R_sms
(mol yr-1),VT_ReactTarget
,description
="total or component R source - sink"R_conc_sms
(mol m-3 yr-1),VT_ReactContributor
,VF_Deriv
,description
="total or component R_conc source - sink"volume
(m3),VT_ReactDependency
,description
="cell volume (as used by total variable)"
do_implicitreservoir_total
R_calc
(mol),VT_ReactTarget
,description
="contributions to total R_calc_conc (NB: a total, not concentration, to generalize to multiphase eqb)"volume
(m3),VT_ReactDependency
,description
="cell volume (as used by total variable)"R
(mol),VT_ReactProperty
,description
="total or component R"R_conc
(mol m-3),VT_ReactContributor
,VF_Total
,description
="total or component R_conc"
PALEOaqchem.ReservoirsAq.ReactionAqConcSum
— TypeReactionAqConcSum
A sum of concentration variables (eg to get an element total)
Parameters
vars_to_add[Vector{String}]
=["2\myvar", "myothervar", "-1\mythirdvar"],default_value
=["2\myvar", "myothervar", "-1\mythirdvar"],description
="vector of variable names to add, eg [2*myvar, myothervar, -1*mythirdvar]"add_to_sum_volume[Bool]
=false,default_value
=false,description
="true to also add to a 'sum' Variable += 'sum_conc * volume"define_sum_volume[Bool]
=false,default_value
=false,description
="only if 'add_to_sum_volume == true': true to also define the 'sum' Variable"
Methods and Variables for default Parameters
do_aqconcsum
sum_conc
(mol m-3),VT_ReactProperty
,description
="sum of specified variables"myvar
(),VT_ReactDependency
,description
=""myothervar
(),VT_ReactDependency
,description
=""mythirdvar
(),VT_ReactDependency
,description
=""
Equilibrium reactions
PALEOaqchem.GenericReactions.ReactionAqEqb
— TypeReactionAqEqb
Define a new equilibrium species N
N + a A <--> b B + c C
[N] = K_eqb'^K_power ([B]^b [C]^c) / ([A]^a)
where to convert density units for K_eqb
:
K_eqb' = K_eqb * rho_ref^K_density_power
The first name in the Reactants
list is the new species concentration N
: other species concentrations in Reactants
and Products
lists must already be defined elsewhere in the model configuration.
The contribution of the new species to element totals or components is defined by the N_components
vector, which may be empty eg to just calculate an Omega or a gas partial pressure etc.
Examples
Gas partial pressure from concentration
solubility_H2:
class: ReactionAqEqb
parameters:
Reactants: ["pH2"]
Products: ["H2_conc"]
K_eqb: 7.8e-1 # mol m-3 atm-1 at 298.15 K (Henry's law coefficent)
K_power: -1.0 # pH2 = H2_conc * K_eqb^-1
variable_attributes:
pH2%units: atm
Compilation of Henry's law coefficients: https://www.henrys-law.org/ which is R. Sander: Compilation of Henry's law constants (version 5.0.0) for water as solvent, Atmos. Chem. Phys., 23, 10901-12440 (2023), doi:10.5194/acp-23-10901-2023
Unit conversions: 1 mol m^-3 Pa-1 = 1.01325e5 mol m-3 atm-1
Parameters
Reactants[Vector{String}]
=["N\conc", "A\conc^2"],default_value
=["N\conc", "A\conc^2"],description
="concentrations or activities of new species followed by other reactants, write powers as X^2 etc"Products[Vector{String}]
=["B\conc", "C\conc"],default_value
=["B\conc", "C\conc"],description
="concentrations or activities of products, write powers as X^2 etc"K_eqb[Float64]
=0.00018629779999999998,default_value
=0.00018629779999999998,description
="equilibrium constant"K_density_power[Float64]
=0.0,default_value
=0.0,description
="multiple K_eqb * rho_ref^K_density_power to convert units: 0.0 for K_eqb in mol m-3, 1.0 for K_eqb in mol kg-1, etc"K_power[Float64]
=-1.0,default_value
=-1.0,description
="exponent of K_eqb"N_components[Vector{String}]
=String[],default_value
=String[],description
="contribution of new species to element or component total eg '["2*TN_calc"]' to add2\*N\_conc\*volume
toTN\_calc
(or empty vector to just define an Omega)"
Methods and Variables for default Parameters
do_aqeqb
N_conc
(mol m-3),VT_ReactProperty
,description
="aqueous concentration or activity"A_conc
(mol m-3),VT_ReactDependency
,description
="aqueous concentration or activity"B_conc
(mol m-3),VT_ReactDependency
,description
="aqueous concentration or activity"C_conc
(mol m-3),VT_ReactDependency
,description
="aqueous concentration or activity"
Kinetic reactions
PALEOaqchem.GenericReactions.ReactionAqKinetic
— TypeReactionAqKinetic
Define a kinetic reaction with rate dependent on concentrations
a A + b B --> c C + d D
Rate (default) is:
R = K * [A] * [B]
where this can be modified to different functional form by defining a vector of Rate_functions
to apply to each concentration variable.
Parameters Reactants
and Products
should be the vectors of stoichiometry * <name> of (total) species to accumulate fluxes into <name>_sms
variables.
Parameter Reactant_concs
should be an empty vector to take default concentration variable names from Reactants
, or a Vector of names to specify concentration species names explicitly (required when eg where Reactants
refers to totals which are partioned into multiple species).
Parameters
Reactants[Vector{String}]
=["A", "2\B"],default_value
=["A", "2\B"],description
="reactant species"Products[Vector{String}]
=["2\C", "D"],default_value
=["2\C", "D"],description
="product species"Reactant_concs[Vector{String}]
=String[],default_value
=String[],description
="names of concentration variables to calculate rate eg '[`"A_conc"]' etc, empty vector to used defaults from 'Reactants' eg 'A_conc', 'B_conc' ..."Rate_functions[Vector{String}]
=String[],default_value
=String[],allowed_values
=["linear", "sqrt", "monod"],description
="functional form for rate function of each concentration (empty vector for default 'linear')"K[Float64]
=NaN,default_value
=NaN,description
="rate constant"K_lim[Float64]
=NaN (mol m-3),default_value
=NaN,description
="limiting concentration for 'monod' rate function"
Methods and Variables for default Parameters
do_aqkinetic
redox_A_B_C_D
(mol yr-1),VT_ReactProperty
,description
="rate variable"A_conc
(mol m-3),VT_ReactDependency
,description
="aqueous concentration or activity"B_conc
(mol m-3),VT_ReactDependency
,description
="aqueous concentration or activity"volume
(m3),VT_ReactDependency
,description
="cell solute volume"
RateStoich_redox_A_B_C_D
redox_A_B_C_D
(mol yr-1),VT_ReactDependency
,description
="rate variable"- [
A_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=redox_A_B_C_D" - [
B_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=redox_A_B_C_D" - [
C_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=redox_A_B_C_D" - [
D_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=redox_A_B_C_D"
totals
redox_A_B_C_D
(mol yr-1),VT_ReactDependency
,description
="rate variable"redox_A_B_C_D_total
(mol yr-1),VT_ReactProperty
,description
="total rate variable"
Precipitation-dissolution reactions
PALEOaqchem.GenericReactions.ReactionAqPrecipDissol
— TypeReactionAqPrecipDissol
Define a precipitation and dissolution reaction for solid S
a A + b B <--> s S + d D
Rate for the precipitation and dissolution reactions are:
R_precip = K_precip * (Ω - 1) (Ω > 1)
R_dissol = K_dissol * S_conc * (1 - Ω) (Ω < 1)
Parameters Reactants
and Products
should be the vectors of stoichiometry * name of (total) species to accumulate fluxes into _sms
variables.
Solid_conc
should be the name of the concentration variable for S, or an empty string to use default 'S_conc'.
Parameters
Reactants[Vector{String}]
=["A", "2\B"],default_value
=["A", "2\B"],description
="reactant species"Products[Vector{String}]
=["S", "0.5\D"],default_value
=["S", "0.5\D"],description
="product species (solid S first)"Solid_conc[String]
="",default_value
="",description
="name of solid S concentration variable (empty string to default to 'S_conc')"K_precip[Float64]
=0.0 (mol m-3 yr-1),default_value
=0.0,description
="rate constant for precipitation reaction"K_dissol[Float64]
=0.0 (yr-1),default_value
=0.0,description
="rate constant for dissolution reaction"dissol_rolloff_conc[Float64]
=0.0 (mol m-3),default_value
=0.0,description
="limiting concentration below which dissolution rate is rolled off to zero as Solid_conc^2"
Methods and Variables for default Parameters
do_aqprecipdissol
precipdissol_A_B_S_D
(mol yr-1),VT_ReactProperty
,description
="rate variable"S_conc
(mol m-3),VT_ReactDependency
,description
="solid concentration or activity"Omega
(),VT_ReactDependency
,description
="saturation state"volume
(m3),VT_ReactDependency
,description
="cell solid phase volume"
RateStoich_precipdissol_A_B_S_D
precipdissol_A_B_S_D
(mol yr-1),VT_ReactDependency
,description
="rate variable"- [
A_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=precipdissol_A_B_S_D" - [
B_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=precipdissol_A_B_S_D" - [
S_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=precipdissol_A_B_S_D" - [
D_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=precipdissol_A_B_S_D"
totals
precipdissol_A_B_S_D
(mol yr-1),VT_ReactDependency
,description
="rate variable"precipdissol_A_B_S_D_total
(mol yr-1),VT_ReactProperty
,description
="total rate variable"
Particulate fluxes
PALEOaqchem.Particle.ReactionParticleDecay
— TypeReactionParticleDecay
Decay (eg organic matter to remineralization) at rate 1.0/decay_timescale
of eg an organic matter dissolved/particulate phase in Reservoir Particle
, to decayflux
. Particle
may have isotope_type
. The Reservoir Particle
may have the :vsink
attribute set to represent a marine sinking particulate phase.
Parameters
decay_timescale[Float64]
=0.5 (yr),default_value
=0.5,description
="particle decay timescale"decay_threshold[Float64]
=-Inf (mol m-3),default_value
=-Inf,description
="particle decay concentration below which decay stops"show_decayrate[Bool]
=false,default_value
=false,description
="true to create a 'decayrate' variable to record decay rate"field_data[DataType]
=PALEOboxes.ScalarData,default_value
=PALEOboxes.ScalarData,allowed_values
=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear],description
="disable / enable isotopes and specify isotope type"
Methods and Variables
do_particle_decay
Particle
(mol),VT_ReactDependency
,description
="Particle amount"Particle_sms
(mol yr-1),VT_ReactContributor
,description
="Particle source-sink"decayflux
(mol yr-1),VT_ReactContributor
,description
="Particle decay flux"
PALEOaqchem.Particle.ReactionFluxToComponents
— TypeReactionFluxToComponents
Distribute a single input_flux
(eg an organic matter flux) to output_fluxnames
components with stoichiometry output_fluxstoich
. input_flux
may have an isotope type (set by field_data
) in which case the isotopic composition will be sent to (usually one) output_fluxname
with ::Isotope
suffix.
Parameters
outputflux_prefix[String]
="",default_value
="",description
="Prefix for output flux names"outputflux_names[Vector{String}]
=["Corg", "N", "P"],default_value
=["Corg", "N", "P"],description
="Suffixes for output flux names. Use ::Isotope suffix to identify a flux with 'isotope_type'"outputflux_stoich[Vector{Float64}]
=[106.0, 16.0, 1.0],default_value
=[106.0, 16.0, 1.0],description
="Stoichiometry for output fluxes relative to input flux"allow_unlinked_outputflux[Bool]
=false,default_value
=false,description
="true to allow (and ignore) unlinked outputflux Variables, false to error if any outputflux Variable is not linked"field_data[DataType]
=PALEOboxes.ScalarData,default_value
=PALEOboxes.ScalarData,allowed_values
=Type[PALEOboxes.ScalarData, PALEOboxes.IsotopeLinear],description
="disable / enable isotopes and specify isotope type"
Methods and Variables for default Parameters
do_flux_to_components
inputflux
(mol yr-1),VT_ReactTarget
,description
="input flux"Corg
(mol yr-1),VT_ReactContributor
,description
="flux Corg"N
(mol yr-1),VT_ReactContributor
,description
="flux N"P
(mol yr-1),VT_ReactContributor
,description
="flux P"
Co-precipitation
PALEOaqchem.CoPrecip.ReactionPACoPrecip
— TypeReactionPACoPrecip
Co-precipitation of P (eg iron-bound phosphorus) with A (eg Fe oxide) formation
P -> P=A
at a rate gamma * A_formation_rate_<n>
, with limitation at low P concentration
P_components
defines solution totals or components that should be modified when P is consumed (eg ["-1*P"]
to remove "P" from solution)
Parameters
A_rate_stoich_factors[Vector{Float64}]
=[1.0],default_value
=[1.0],description
="stoichiometry factor to multiply each A formation rate variable to convert to mol A"gamma[Float64]
=0.15 (mol/mol),default_value
=0.15,description
="P:A molar ratio"P_limit[Float64]
=1.0e-6 (mol m-3),default_value
=1.0e-6,description
="limiting P concentration below which co-precipitation is inhibited"P_components[Vector{String}]
=["-1\P", "TAlk"],default_value
=["-1\P", "TAlk"],description
="totals or components that should be modified when P is consumed from solution"
Methods and Variables for default Parameters
do_PA_coprecip
rate_PA_coprecip
(mol P yr-1),VT_ReactProperty
,description
="rate of P co-precipitation"P_conc
(mol m-3),VT_ReactDependency
,description
="P concentration"A_formation_rate_1
(mol m-3 yr-1),VT_ReactDependency
,description
="substance A formation rate"
RateStoich_rate_PA_coprecip
rate_PA_coprecip
(mol P yr-1),VT_ReactDependency
,description
="rate of P co-precipitation"- [
P_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_coprecip" - [
TAlk_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_coprecip" - [
PA_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_coprecip"
totals
rate_PA_coprecip
(mol P yr-1),VT_ReactDependency
,description
="rate of P co-precipitation"rate_PA_coprecip_total
(mol P yr-1),VT_ReactProperty
,description
="total rate of P co-precipitation"
PALEOaqchem.CoPrecip.ReactionPARelease
— TypeReactionPARelease
Release of P (eg iron-bound phosphorus) with A (eg Fe oxide) destruction
P=A -> P
at a rate Prelease = theta * A_destruction_rate_<n>
, where theta = PA_conc / A_conc
P_components
defines totals or components that should be modified when P is released (eg ["P"]
to return P
to solution).
Parameters
A_rate_stoich_factors[Vector{Float64}]
=[1.0],default_value
=[1.0],description
="stoichiometry factor to multiply each A destruction rate variable to convert to mol A"P_components[Vector{String}]
=["P", "-1\TAlk"],default_value
=["P", "-1\TAlk"],description
="totals or components that should be modified when P is released"
Methods and Variables for default Parameters
do_PA_release
rate_PA_release
(mol P yr-1),VT_ReactProperty
,description
="rate of coprecipitated P dissolution"PA_conc
(mol m-3),VT_ReactDependency
,description
="adsorbed P concentration"A_conc
(mol m-3),VT_ReactDependency
,description
="adsorbant concentration"PA_theta
(mol/mol),VT_ReactProperty
,description
="P:A molar ratio of adsorbed of coprecipitated P"A_destruction_rate_1
(mol m-3 yr-1),VT_ReactDependency
,description
="substance A destruction rate"
RateStoich_rate_PA_release
rate_PA_release
(mol P yr-1),VT_ReactDependency
,description
="rate of coprecipitated P dissolution"- [
PA_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_release" - [
P_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_release" - [
TAlk_sms
] (mol yr-1),VT_ReactContributor
,description
="generated by RateStoich rate=rate_PA_release"
totals
rate_PA_release
(mol P yr-1),VT_ReactDependency
,description
="rate of coprecipitated P dissolution"rate_PA_release_total
(mol P yr-1),VT_ReactProperty
,description
="total rate of coprecipitated P dissolution"